ex12
Lindsey Pekurny and Reed Scott
4/19/2022
Question 1
- Build ensemble map products of:
- SDHM mean and sd probabilities
- SDHM concordance maps of 5, 3, and union of all SDHM overlaps
- “Clip” all these maps by the bounding boxes created earlier (see Module 2.3.3 for refresher, if needed)
- Output these map products as
.imgfiles
#path.root <- "D:/OneDrive - University of Vermont/Classes/Spring2022/sdhmR/sdhmR-V2022.1" #Reed Laptop Path
path.root <- "C:/Users/14842/Documents/SDHM/sdhmR-V2022.1" #Lindsey Path
path.mod <- paste(path.root, "/data/exercise/traindat", sep = "")
path.preds <- paste(path.root, '/data/exercise/preds', sep = '')
path.maps <- paste(path.root, '/data/exercise/maps', sep = '')
path.gis <- paste(path.root, "/data/gis_layers", sep = "")
path.figs <- paste(path.root, "/powerpoints/figures", sep = "")load libraries
library(raster)
options("rgdal_show_exportToProj4_warnings" = "none") # run this before library(rgdal)
library(rgdal)
library(gam)
library(randomForest)
library(dismo)
library(gbm)
library(sf)START INITIALIZATION OF DATA STRUCTURES
load SDM models. I renamed the models and the threshold cuts to make it more clear raster stack of predictors assumes rasters all have same projection, extent & resolution
setwd(path.preds)
pers.list <- list.files(pattern = ".img$") # list of .img files; $ strips extra
pers.list # examine## [1] "etpt_5.img" "etpt_6.img" "etpt_sprin.img" "exp1nrm.img"
## [5] "exp3nrm.img" "exp5nrm.img" "mind_yr_av.img" "prad_sw_di.img"
## [9] "prec_w_hal.img" "prec_winte.img" "rough_1k.img" "tave_s_hal.img"
## [13] "tave_sprin.img" "tmax_s_hal.img" "tmax_summe.img" "topos.img" pers.dom <- stack(pers.list) # build raster stack
pers.dom # examine stack## class : RasterStack
## dimensions : 1518, 1478, 2243604, 16 (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -114.6167, -102.3, 29.49167, 42.14167 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : etpt_5, etpt_6, etpt_sprin, exp1nrm, exp3nrm, exp5nrm, mind_yr_av, prad_sw_di, prec_w_hal, prec_winte, rough_1k, tave_s_hal, tave_sprin, tmax_s_hal, tmax_summe, ...
## min values : 0.000000, 69.000000, 0.000000, -1152.000000, -812.000000, -715.000000, -2294.333252, 6411.166016, 3.666667, 3.333333, 0.000000, 16.166666, -86.000000, 63.833332, 100.000000, ...
## max values : 100.00000, 104.00000, 90.33334, 1246.00000, 945.00000, 895.00000, 249.33333, 14441.50000, 136.83333, 140.33333, 243.85391, 308.33334, 226.66667, 392.33334, 424.00000, ... names(pers.dom) ## [1] "etpt_5" "etpt_6" "etpt_sprin" "exp1nrm" "exp3nrm"
## [6] "exp5nrm" "mind_yr_av" "prad_sw_di" "prec_w_hal" "prec_winte"
## [11] "rough_1k" "tave_s_hal" "tave_sprin" "tmax_s_hal" "tmax_summe"
## [16] "topos"# load SDM models. I renamed the models and the threshold cuts to make it more clear
setwd(paste(path.mod, sep = ""))
load('ex7.RData')
LR.cut <- mod.cut$glm.cvpred ### Assuming this is the cut point? In his example it's just one value
LR <- mod2.LR
load('ex8.RData')
GAM.cut <- cut$pred
GAM <- mod1.GAM
load('ex9.RData')
MAX <- mod1.MAX
MAX.cut <- mod.cut$spec_sens
load('ex10.RData')
RF <- pers.RF
RF.cut <- mod.cut$persRF.pred
load('ex11.RData')
BRT <- pers.BRT
BRT.cut <- mod.cut$pred
### Building a list of cut points.
cut.list <- c(LR.cut,
GAM.cut,
MAX.cut,
RF.cut,
BRT.cut
)
cut.list # threshold cuts?## [1] 0.5077240 0.6200000 0.4894119 0.3400000 0.5077240# raster stack of predictors assumes rasters all have same projection, extent & resolution
setwd(path.preds)
pers.list <- list.files(pattern = ".img$") # list of .img files; $ strips extra
pers.list # examine## [1] "etpt_5.img" "etpt_6.img" "etpt_sprin.img" "exp1nrm.img"
## [5] "exp3nrm.img" "exp5nrm.img" "mind_yr_av.img" "prad_sw_di.img"
## [9] "prec_w_hal.img" "prec_winte.img" "rough_1k.img" "tave_s_hal.img"
## [13] "tave_sprin.img" "tmax_s_hal.img" "tmax_summe.img" "topos.img" pers.dom <- stack(pers.list) # build raster stack
pers.dom # examine stack## class : RasterStack
## dimensions : 1518, 1478, 2243604, 16 (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -114.6167, -102.3, 29.49167, 42.14167 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : etpt_5, etpt_6, etpt_sprin, exp1nrm, exp3nrm, exp5nrm, mind_yr_av, prad_sw_di, prec_w_hal, prec_winte, rough_1k, tave_s_hal, tave_sprin, tmax_s_hal, tmax_summe, ...
## min values : 0.000000, 69.000000, 0.000000, -1152.000000, -812.000000, -715.000000, -2294.333252, 6411.166016, 3.666667, 3.333333, 0.000000, 16.166666, -86.000000, 63.833332, 100.000000, ...
## max values : 100.00000, 104.00000, 90.33334, 1246.00000, 945.00000, 895.00000, 249.33333, 14441.50000, 136.83333, 140.33333, 243.85391, 308.33334, 226.66667, 392.33334, 424.00000, ... names(pers.dom)## [1] "etpt_5" "etpt_6" "etpt_sprin" "exp1nrm" "exp3nrm"
## [6] "exp5nrm" "mind_yr_av" "prad_sw_di" "prec_w_hal" "prec_winte"
## [11] "rough_1k" "tave_s_hal" "tave_sprin" "tmax_s_hal" "tmax_summe"
## [16] "topos"# LR prediction and classification
setwd(paste(path.maps,
sep = ""))
names(pers.dom) <- paste(names(pers.dom),'img')
modFprob.LR <- predict(pers.dom,
LR,
filename = "modFprob.GLM.img",
type = "response",
fun = predict,
index = 2,
overwrite = T)
modFclass.LR <- reclassify(pers.prob,
filename = "modFclass.LR.img",
c(0,
mod.cut[[2]],
0,
mod.cut[[2]],
1,
1),
overwrite=TRUE)# giggle plots
par(mfrow = c(1, 2))
plot(modFprob.LR,
axes = T,
main = "LR probability map") # plot probability map
plot(modFclass.LR,
axes = T,
main = "LR classfied map") # plot classified map
par(mfrow = c(1, 1))# GAM prediction and classification
library(gam)
modFprob.GAM <- predict(pers.dom,
GAM,
filename = "modFprob.GAM.img",
type = "response",
fun = predict,
index = 2,
overwrite = T)
modFclas.GAM <- reclassify(modFprob.GAM,
filename = "modFclas.GAM.img",
(c(0,
GAM.cut,
0,
GAM.cut,
1,
1)),
overwrite = T) # giggle plots
par(mfrow = c(1, 2))
plot(modFprob.GAM,
axes = T,
main = "GAM probability map") # plot probability map
plot(modFclas.GAM,
axes = T,
main = "GAM classfied map") # plot classified map
par(mfrow = c(1, 1)) # RF prediction and classification
library(randomForest) # load library
modFprob.RF <- predict(pers.dom,
RF,
filename = "modFprob.RF.img",
type = "prob",
fun = predict,
index = 2,
overwrite = T) # prob map
modFclas.RF <- reclassify(modFprob.RF,
filename = "modFclas.RF.img",
(c(0,
RF.cut,
0,
RF.cut,
1,
1)),
overwrite = T) # class map # giggle maps
par(mfrow = c(1, 2))
plot(modFprob.RF,
axes = T,
main = "RF probability map") # plot probability map
plot(modFclas.RF,
axes = T,
main = "RF classfied map") # plot classified map
par(mfrow = c(1, 1)) library(dismo)
library(gbm)
modFprob.BRT <- predict(pers.dom,
BRT,
n.trees = BRT$gbm.call$best.trees,
type = "response",
filename = "modFprob.BRT.img",
overwrite = T) # prob map
modFclas.BRT <- reclassify(modFprob.BRT,
c(0, BRT.cut,
0, BRT.cut,
1,
1),
filename = "modFclas.BRT.img",
overwrite = T) # clas map # giggle plots
par(mfrow = c(1, 2))
plot(modFprob.BRT,
axes = T,
main = "BRT probability map") # plot probability map
plot(modFclas.BRT, axes = T,
main = "BRT classfied map") # plot classified map
par(mfrow = c(1, 2))
# save BRT plots if desired
setwd(path.figs)
savePlot(filename = "mod06fig04.pdf",
type = "pdf")## Error in savePlot(filename = "mod06fig04.pdf", type = "pdf"): can only copy from 'windows' devices ## MAXENT prediction and classification
setwd(path.maps)
library(dismo) # load library
modFprob.MAX <- predict(MAX,
pers.dom,
filename = "modFprob.MAX.img",
overwrite = T)
modFclas.MAX <- reclassify(modFprob.MAX,
filename = "modFclas.MAX.img",
c(0,
MAX.cut,
0,
MAX.cut,
1,
1),
overwrite = T) # clas map
# giggle plots
par(mfrow = c(1, 2))
plot(modFprob.MAX, axes = T, main = "MAX probability map") # plot probability map
plot(modFclas.MAX, axes = T, main = "MAX classfied map") # plot classified map
par(mfrow = c(1, 1))
# save MAX plots if desired
setwd(path.figs)
savePlot(filename = "mod06fig05.pdf", type = "pdf")## Error in savePlot(filename = "mod06fig05.pdf", type = "pdf"): can only copy from 'windows' devicesEND PREDICTION PROBABILITY AND CLASSIFED MAPS START ENSEMBLE PROCESS # load probability & classified maps; create stacks
setwd(paste(path.maps, sep = ""))
pr.list <- unlist(unique(strsplit(list.files(pattern = "Fprob"),
".aux.xml")))
pr.list # examine## [1] "modFprob.BRT.img" "modFprob.GAM.img" "modFprob.GLM.img" "modFprob.MAX.img"
## [5] "modFprob.RF.img" prob.dom <- stack(pr.list) # raster stack prob maps
prob.dom # examine## class : RasterStack
## dimensions : 1518, 1478, 2243604, 5 (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -114.6167, -102.3, 29.49167, 42.14167 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : modFprob.BRT, modFprob.GAM, modFprob.GLM, modFprob.MAX, modFprob.RF
## min values : 2.186834e-01, 2.220446e-16, 2.220446e-16, 1.048621e-05, 0.000000e+00
## max values : 0.6525138, 0.9871084, 0.9911523, 0.7762663, 1.0000000 # cl.list <- list.files(pattern = "Fclas")
cl.list <- unlist(unique(strsplit(list.files(pattern = "Fclas"),
".aux.xml")))
cl.list # examine## [1] "modFclas.BRT.img" "modFclas.GAM.img" "modFclas.LR.img" "modFclas.MAX.img"
## [5] "modFclas.RF.img" "modFclass.LR.img" clas.dom <- stack(cl.list) # raster stack classified maps
clas.dom # examine## class : RasterStack
## dimensions : 1518, 1478, 2243604, 6 (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -114.6167, -102.3, 29.49167, 42.14167 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : modFclas.BRT, modFclas.GAM, modFclas.LR, modFclas.MAX, modFclas.RF, modFclass.LR
## min values : 0, 0, 0, 0, 0, 0
## max values : 1, 1, 1, 1, 1, 1 # standardize all prediction maps 0-1.0
maxValue(prob.dom) # extract max value for each prob map## [1] 0.6525138 0.9871084 0.9911523 0.7762663 1.0000000 layers <- {} # initialize (empty) list of raster layers
for (i in 1:length(names(prob.dom))) {
m1 <- prob.dom[[i]] # get a prob map
m2 <- 1/maxValue(prob.dom[[i]]) * prob.dom[[i]] # standardize all probs to max=1
m3 <- unlist(strsplit(names(prob.dom[[i]]), "[.]")) # split prob layer name apart
names(m2) <- paste(m3[1], "STD.", m3[2], sep = "") # assign name to raster value
assign(paste(m3[1], "STD.", m3[2], sep = ""), m2) # assign new name to standardized layer
layers <- c(layers, get(paste(m3[1], "STD.", m3[2], sep = "")))
}
probSTD.dom <- stack(layers)
maxValue(probSTD.dom) # extract max value for each prob map; standardized ?## [1] 1 1 1 1 1 # save stacks
setwd(path.mod)
save(prob.dom,
probSTD.dom,
clas.dom,
file = "ensemble.dom.RData") ######################################################
# descriptive stats on raster maps: mean & sum
prob.mean <- mean(prob.dom) # mean prob map
# access min max values of descriptive stat raster
maxValue(prob.mean) # max pr.mean## [1] 0.8468679 minValue(prob.mean) # min pr.mean## [1] 0.04380246 probSTD.mean <- mean(probSTD.dom) # mean standardized prob map
clas.sum <- sum(clas.dom) # sum of models by cell
# giggle plots: mean & sum
par(mfrow = c(1, 2))
plot(prob.mean,
axes = T,
main = "MEAN probability map") # plot probability map
plot(st_geometry(pegr6.mahog),
add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame),
add = T) # make pretty## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found plot(clas.sum,
axes = T,
main = "Concordance CLASS map: Ramp") # plot concordance map
plot(st_geometry(pegr6.mahog),
add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame),
add = T) # make pretty ## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found par(mfrow = c(1, 1))
# save plots if desired
setwd(path.figs)
savePlot(filename = "mod06fig06.pdf", type = "pdf")## Error in savePlot(filename = "mod06fig06.pdf", type = "pdf"): can only copy from 'windows' devices ####
# arithmetic on raster maps: multiplication & subsetting
prob.mean1 <- (clas.sum > 0) * prob.mean # mean where class=1
clas.sum1 <- clas.sum > 0 # sum where No. models >=1
# giggle plots: means X classifications
par(mfrow = c(1, 2))
plot(prob.mean1,
axes = T,
main = "MEAN probability map: Presence=1") # plot probability map
plot(st_geometry(pegr6.mahog),
add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame),
add = T) # make pretty## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found plot(clas.sum1, axes = T,
main = "Concordance CLASS map: Union") # plot concordance map
plot(st_geometry(pegr6.mahog), add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame), add = T) # make pretty ## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found par(mfrow = c(1, 1))
####
# calculate measures of variation
prob.sd <- calc(prob.dom, sd) # sd prob map
prob.var <- calc(prob.dom, var) # var prob map
prob.cv <- (prob.sd/prob.mean)*100 # coefficient of variation
# giggle plots: sd and var
par(mfrow = c(1, 2))
plot(prob.sd,
axes = T,
main = "SD probability map") # plot classified map
plot(st_geometry(pegr6.mahog), add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame), add = T) # make pretty## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found plot(prob.var, axes = T, main = "VAR probability map") # plot classified map
plot(st_geometry(pegr6.mahog), add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame), add = T) # make pretty ## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found# make pretty
# save plots if desired #setwd(path.figs) #savePlot(filename = “mod06fig13.pdf”, type = “pdf”) ######## END ENSEMBLE INTERPRETATION ################################################################################ ########################################################################### ## Question #2 # Calculate the frequencies of “presence” in each of the 5 SDHMs
#### # some summary statistics: frequencies
clas.freqM <- freq(clas.dom) # [0,1] freqs by clas map; rtns list
clas.freqM[c(1,5)] # examine freqs for 1st 2 models## $modFclas.BRT
## value count
## [1,] 0 1326676
## [2,] 1 865593
## [3,] NA 51335
##
## $modFclas.RF
## value count
## [1,] 0 1409045
## [2,] 1 783224
## [3,] NA 51335 clas.freqS <- data.frame(freq(clas.sum)) # freqs of concordance of models
clas.freqS # examine; value is concordance freq## value count
## 1 0 1240510
## 2 1 74298
## 3 2 111049
## 4 3 115015
## 5 4 39812
## 6 5 144216
## 7 6 467369
## 8 NA 51335 clas.freqS$count[6]/sum(clas.freqS$count[2:6]) # prop. models concordance=5 ## [1] 0.297727######## START ENSEMBLE INTERPRETATION #### # 6 prob map plots
par(mfrow = c(2, 3))
plot(modFprob.LR,
axes= T,
main = "LR probability map") # LR prob map
plot(st_geometry(pegr6.mahog),
add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame),
add = T) # make pretty## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found plot(modFprob.GAM,
axes = T,
main = "GAM probability map") # GAM prob map
plot(st_geometry(pegr6.mahog), add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame), add = T) # make pretty## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found plot(modFprob.MAX,
axes = T,
main = "MAX probability map") # MAX prob map
plot(st_geometry(pegr6.mahog), add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame), add = T) # make pretty## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found plot(modFprob.RF,
axes = T,
main = "RF probability map") # RF prob map
plot(st_geometry(pegr6.mahog), add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame), add = T) # make pretty## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found plot(modFprob.BRT,
axes = T,
main = "BRT probability map") # BRT prob map
plot(st_geometry(pegr6.mahog), add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame), add = T) # make pretty ## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found plot(prob.mean1,
axes = T,
main = "MEAN probability map: \nPresence=1") # MEAN prob map
plot(st_geometry(pegr6.mahog), add = T)## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found plot(st_geometry(pegr6.frame), add = T) # make pretty ## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found#### # 6 classified maps
par(mfrow = c(2, 3))
plot(modFclass.LR,
axes = T,
main = "LR classfied map") # LR classified map
plot(modFclas.GAM, axes = T, main = "GAM classfied map") # GAM classified map
plot(modFclas.MAX, axes = T, main = "MAX classfied map") # MAX classified map
plot(modFclas.RF, axes = T, main = "RF classfied map") # RF classified map
plot(modFclas.BRT, axes = T, main = "BRT classfied map") # BRT classified map
plot(clas.sum, axes = T, main = "Concordance CLASS \nmap: Ramp") # plot concordance map
# save plots if desired
#setwd(path.figs)
#savePlot(filename = "mod06fig11.pdf", type = "pdf") ####
# 6 classification union maps
clas.sumU <- clas.sum > 0 # sum where No. models =union
clas.sum1 <- clas.sum > 1 # sum where No. models =2+
clas.sum2 <- clas.sum > 2 # sum where No. models =3+
clas.sum3 <- clas.sum > 3 # sum where No. models =4+
clas.sum4 <- clas.sum > 4 # sum where No. models =5
# giggle plots
par(mfrow = c(2, 3))
plot(clas.sum, axes = T, main = "Concordance CLASS \nmap: Ramp") # RAMP concordance map
plot(clas.sumU, axes = T, main = "Concordance CLASS \nmap: Union") # UNION concordance map
plot(clas.sum1, axes = T, main = "Concordance CLASS \nmap: 2+") # 2+ concordance map
plot(clas.sum2, axes = T, main = "Concordance CLASS \nmap: 3+") # 3+ concordance map
plot(clas.sum3, axes = T, main = "Concordance CLASS \nmap: 4+") # 4+ concordance map
plot(clas.sum4, axes = T, main = "Concordance CLASS \nmap: 5") # 5 concordance map
# save plots if desired #setwd(path.figs) #savePlot(filename = “mod06fig12.pdf”, type = “pdf”)
####
# some variance analyses
q1 <- summary(prob.sd)[2] # 1st quartile
q2 <- summary(prob.sd)[3] # 2nd quartile
q3 <- summary(prob.sd)[4] # 3rd quartile
prob.sd1q <- (prob.sd >= q1) * clas.sum1 # where 2+ models agree
prob.sd2q <- (prob.sd >= q2) * clas.sum1 # where 2+ models agree
prob.sd3q <- (prob.sd >= q3) * clas.sum1 # where 2+ models agree
# giggle plots
par(mfrow = c(2, 2))
plot(prob.sd, axes = T, main = "SD probability map") # SD prob map
plot(st_geometry(pegr6.mahog), add = T)
plot(st_geometry(pegr6.frame), add = T) # make pretty
plot(prob.sd1q, axes = Tmain = "1st quartile SD map") # 1stQ prob map
plot(st_geometry(pegr6.mahog), add = T)
plot(st_geometry(pegr6.frame), add = T) # make pretty
plot(prob.sd2q, axes = T, main = "2nd quartile SD map") # 2ndQ prob map
plot(st_geometry(pegr6.mahog), add = T)
plot(st_geometry(pegr6.frame), add = T) # make pretty
plot(prob.sd3q, axes = T, main = "3rd quartile SD map") # 3rdQ prob map
plot(st_geometry(pegr6.mahog), add = T)
plot(st_geometry(pegr6.frame), add = T)## Error: <text>:16:32: unexpected '='
## 15:
## 16: plot(prob.sd1q, axes = Tmain =
## ^Question #3
#Tally the frequencies of concordance after “clipping” and compare with above
#### clipping...
# Need to clip pers.DOM
clas.freqS <- data.frame(freq(clas.sum)) # freqs of concordance of models
clas.freqS # examine; value is concordance freq## value count
## 1 0 1240510
## 2 1 74298
## 3 2 111049
## 4 3 115015
## 5 4 39812
## 6 5 144216
## 7 6 467369
## 8 NA 51335clas.freqS$count[6]/sum(clas.freqS$count[2:6])## [1] 0.297727 # 6 classification union maps
clas.sumU.freq <- data.frame(freq(clas.sumU)) # sum where No. models =union
clas.sum1.freq <- data.frame(freq(clas.sum1)) # sum where No. models =2+
clas.sum2.freq <- data.frame(freq(clas.sum2)) # sum where No. models =3+
clas.sum3.freq <- data.frame(freq(clas.sum3)) # sum where No. models =4+
clas.sum4.freq <- data.frame(freq(clas.sum4))# sum where No. models =5
clas.sumU.freq$count## [1] 1240510 951759 51335 clas.sum1.freq$count ## [1] 1314808 877461 51335 clas.sum2.freq$count ## [1] 1425857 766412 51335 clas.sum3.freq$count ## [1] 1540872 651397 51335 clas.sum4.freq$count## [1] 1580684 611585 51335# Concordance frequencies decline as maps are layered upon one another.Question #4
#* Save your data as R objects: # * All ensemble prediction maps as .img format #* Save these R objects in a .RData file as well
setwd(path.ex)## Error in setwd(path.ex): object 'path.ex' not foundsave(pers.BRT,pers.BRT2, file = 'ex12.RData')
save(pers.class.BRT, file = 'persClassBRT.img')## Error in save(pers.class.BRT, file = "persClassBRT.img"): object 'pers.class.BRT' not foundsave(pers.prob.BRT, file = 'persProbBRT.img')## Error in save(pers.prob.BRT, file = "persProbBRT.img"): object 'pers.prob.BRT' not found